home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / GAMMA / GAMMA.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1991-03-17  |  4.6 KB  |  150 lines
  1. PROGRAM GammaTest;
  2. {$N+}
  3. {$E+}
  4. VAR
  5.    Quit:BOOLEAN;
  6.    z:EXTENDED;
  7.  
  8. FUNCTION Gamma(x:EXTENDED):EXTENDED;
  9. { Compute gamma to a fair number of significant digits.  Written by
  10.   Tom Hardy 16 March 1991 }
  11. CONST
  12.    Big=1000.0;
  13.    MaxReal=1.0E100;
  14.  
  15. FUNCTION NearZero(x:EXTENDED):BOOLEAN;
  16. BEGIN { NearZero }
  17. NearZero:=(ABS(x)<1.0E-18)
  18. END; { NearZero }
  19.  
  20. FUNCTION NearInteger(x:EXTENDED):BOOLEAN;
  21. BEGIN { NearInteger }
  22. NearInteger:=NearZero(FRAC(x)) OR NearZero(1-ABS(FRAC(x)))
  23. END; { NearInteger }
  24.  
  25. FUNCTION Gamma2(z:EXTENDED):EXTENDED;
  26. { Taylor expansion for Γ(z+3).  From Mathematical Functions and Their
  27.   Approximations by Yudell L. Luke, Academic Press 1975 }
  28. CONST
  29.    Max=41;
  30.    b:ARRAY [0.. Max] OF EXTENDED=(2.00000000000000000000,
  31.                                   1.84556867019693427879,
  32.                                   1.24646499595134652897,
  33.                                   0.57499416892061222755,
  34.                                   0.23007494075411406302,
  35.                                   0.07371504661602386878,
  36.                                   0.02204110936751696733,
  37.                                   0.00544875407582030942,
  38.                                   0.00135522086023943520,
  39.                                   0.00026478566304549638,
  40.                                   0.00006120306281920073,
  41.                                   0.00000850557917488135,
  42.                                   0.00000240617724013144,
  43.                                   0.00000008802390990648,
  44.                                   0.00000011422276453422,
  45.                                  -0.00000001631475210083,
  46.                                   0.00000000862349738998,
  47.                                  -0.00000000244110402524,
  48.                                   0.00000000087291506387,
  49.                                  -0.00000000028390295131,
  50.                                   0.00000000009560889805,
  51.                                  -0.00000000003178270013,
  52.                                   0.00000000001061093470,
  53.                                  -0.00000000000353684281,
  54.                                   0.00000000000117938189,
  55.                                  -0.00000000000039318677,
  56.                                   0.00000000000013108214,
  57.                                  -0.00000000000004369853,
  58.                                   0.00000000000001456734,
  59.                                  -0.00000000000000485607,
  60.                                   0.00000000000000161876,
  61.                                  -0.00000000000000053961,
  62.                                   0.00000000000000017987,
  63.                                  -0.00000000000000005996,
  64.                                   0.00000000000000001999,
  65.                                  -0.00000000000000000666,
  66.                                   0.00000000000000000222,
  67.                                  -0.00000000000000000074,
  68.                                   0.00000000000000000025,
  69.                                  -0.00000000000000000008,
  70.                                   0.00000000000000000003,
  71.                                  -0.00000000000000000001);
  72. VAR
  73.    n:INTEGER;
  74.    Sum,Powerz:EXTENDED;
  75. BEGIN { Gamma2 }
  76. Sum:=0.0;
  77. Powerz:=1.0;
  78. FOR n:=0 TO Max DO
  79.    BEGIN
  80.    Sum:=Sum+b[n]*Powerz;
  81.    Powerz:=Powerz*z
  82.    END;
  83. Gamma2:=Sum
  84. END; { Gamma2 }
  85.  
  86. FUNCTION Gamma1(z:EXTENDED):EXTENDED;
  87. { If z>=k then map z into the range [1-k,k] }
  88. CONST
  89.    k=1.5;
  90. VAR
  91.    n,i:INTEGER;
  92.    Product:EXTENDED;
  93. BEGIN { Gamma1 }
  94. IF (z>=k) THEN
  95.    BEGIN
  96.    Product:=1.0;
  97.    n:=TRUNC(z-k)+1;
  98.    FOR i:=-n TO -1 DO
  99.       Product:=Product*(i+z);
  100.    Gamma1:=Gamma2(z-n)*Product/((z-n)*(z-n+1)*(z-n+2))
  101.    END
  102. ELSE
  103.    Gamma1:=Gamma2(z)/(z*(1+z)*(2+z))
  104. END; { Gamma1 }
  105.  
  106. BEGIN { Gamma }
  107. IF ((x>Big) OR (x<=-Big)) THEN
  108.    BEGIN
  109.    WRITELN('Invalid gamma function argument');
  110.    Halt(1)
  111.    END
  112. ELSE
  113.    IF (x<=0.0) THEN
  114.       BEGIN
  115.       IF NearInteger(x) THEN
  116.          Gamma:=MaxReal
  117.       ELSE
  118.          Gamma:=Pi/(Gamma1(1-x)*SIN(Pi*x))
  119.       END
  120.    ELSE
  121.       Gamma:=Gamma1(x)
  122. END; { Gamma }
  123.  
  124. PROCEDURE GetReal(Str:STRING;VAR x:EXTENDED;VAR Quit:BOOLEAN);
  125. VAR
  126.    Valid:BOOLEAN;
  127.    Code:INTEGER;
  128.    s:STRING;
  129. BEGIN { GetReal }
  130. REPEAT
  131.    WRITE(Str);
  132.    READLN(s);
  133.    Quit:=(s='');
  134.    Val(s,x,Code);
  135.    Valid:=(Code=0)
  136. UNTIL (Valid OR Quit)
  137. END; { GetReal }
  138.  
  139. BEGIN { GammaTest }
  140. WRITELN('Compute the gamma function to a fair number of significant digits');
  141. WRITELN;
  142. GetReal('Enter z, [ENTER] to quit : ',z,Quit);
  143. WHILE (NOT Quit) DO
  144.    BEGIN
  145.    WRITELN('Γ(z)= ',Gamma(z));
  146.    WRITELN;
  147.    GetReal('Enter z, [ENTER] to quit : ',z,Quit);
  148.    END
  149. END. { GammaTest }
  150.